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We present a position Langevin equation for overdamped particle motion on rough two- 
dimensional surfaces. A Brownian Dynamics algorithm is suggested to evolve this equation nu- 
merically, allowing for the prediction of effective (projected) diffusion coefficients over corrugated 
surfaces. In the case of static surface roughness, we find that a simple area-scaling prediction for the 
projected diffusion coefficient leads to seemingly quantitative agreement with numerical results. To 
study the effect of dynamic surface evolution on the diffusive process, we consider particle diffusion 
over a thermally fiuctuating elastic membrane. Surface fluctuation has the effect of increasing the 
effective diffusivity toward a limiting annealed-surface value discussed previously. We argue that 
protein motion over cell surfaces spans a variety of physical regimes, making it impossible to identify 
a single approximation scheme appropriate to all measurements of interest. 



I. INTRODUCTION 

Diffusive processes are commonplace in many systems 
of chemical, physical and biological interest [H |5J [31 H] , a 
fact commonly attributed to inherent randomness at the 
molecular level of dynamics (thermal or otherwise) and 
the wide range of applicability of the central limit the- 
orem [3]. In the context of biological systems, diffusion 
plays a key role in many of the processes of life at the 
cellular and sub-cellular levels ^ i5j . In particular, the 
diffusion of molecules on the surface of cells has received 
widespread attention for many years [51 [TJ [HI IHl [10] , both 
due to the clear consequences this motion has for cellular 
functioning as well as the availability of suitable exper- 
imental techniques to make quantitative measurements. 
Developments in experimental techniques such as single- 
particle tracking [10 , fluorescence photo-bleaching recov- 
ery [11] and nuclear magnetic resonance (NMR) have 
made it possible to gain direct insight into the structure 
and dynamics of lipid membranes by investigating the 
lateral diffusion of lipids and integral membrane proteins 

Biological membranes are typically not flat, which 
leads to ambiguity in assigning the "diffusion coefficient" 
for a molecule on such a surface. Do we mean the in- 
trinsic curvilinear diffusion coefficient, Z?o, that would 
be expected if the molecule were moving over a truly flat 
surface, or do we mean the effective diffusion coefficient, 
D, reflecting projected motion of the particle onto a ref- 
erence plane (see Fig. [ij? Most experiments, including 
optical techniques and NMR [THl [H] [H], directly mea- 
sure D, while traditional theories provide a route to es- 
timation of Dq [m [15] . The connection between D and 
Dq is the focus of this paper. 

It is well appreciated that confining a particle to a 
two-dimensional (2D) surface within a three-dimensional 
(3D) space has the effect of lowering the effective dif- 
fusion coefficient, D, relative to the curvilinear value. 



Do [Tni[ni[THl[ini[2ni[lIl[12]- Recently, this ge- 
ometric effect has seen renewed interest in the con- 
text of diffusion over lipid bilayers and biomembranes 
[221 [m [551 US]- Reflecting the diversity of membrane 
structures present in biology, some studies have focused 
interest on the case of a particle diffusing over a static 
membrane [THl [231 [2S]) while others have concentrated 
on membranes with fast dynamic fluctuations [531 155] . 
For certain biological structures, the static membrane 
(quenched disorder) limit is probably a good one. For 
example, the endoplasmic reticulum (ER) 5 and the lat- 
eral cortex of auditory outer hair cells ^271, ^28^ are oddly 
shaped in a predominantly static fashion due to topol- 
ogy of the membrane and interactions with other cellular 
components. Other systems, like reconstituted unilamel- 
lar vesicles 29J and the plasma membrane of the red 
blood cell [3.0, are expected to fluctuate rapidly, which 
motivated Reister and Seifert to consider the fast shape- 
fluctuation (annealed disorder) limit in the context of 
membrane protein diffusion 24J. The annealed disorder 
limit is especially attractive for theoretical study, since 
spatial correlations of the membrane surface drop out in 
this limit and analytical results are readily obtained [16]. 
More generally, membrane systems may fall in an inter- 
mediate regime where membrane motion and diffusive 
transport share comparable characteristic time scales. 

While the annealed disorder limit is well understood, 
the opposite limit of quenched disorder is less established. 
One of the studies mentioned above [2B] introduces an 
advanced computational strategy to analyze real experi- 
mental data for static configurations of ER membranes, 
but avoids detailed theoretical interpretation of the re- 
sults. The static surface study by King '^S], is plagued 
by an imprecise treatment of the diffusive process, which 
unfortunately leads to spurious results. Gustafsson and 
Halle [16J suggested a number of plausible approxima- 
tions and rigorous bounds for the behavior of D over a 
static surface, and hypothesized that an effective medium 
approximation (EM A) ^T] would provide the closest cor- 
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respondence with reality, but provided no direct evidence 
to support this conclusion. 

Outside the field of biological physics, general simula- 
tion strategies have been suggested for the evolution of 
diffusive particles over curved surfaces of arbitrary com- 
plexity [IHl HO] . While these methods can be applied to 
membrane systems, in many cases of interest they are 
more elaborate and computationally expensive than nec- 
essary. Commonly, biomcmbrane surfaces are described 
in the Monge gauge. The shape of the membrane surface 
is described by a single- valued function, h{x,y), of lat- 
eral coordinates, x and y h{x,y) provides the local 
height of the membrane above some arbitrary reference 
plane (see Fig. [ij. Given this parameterization, it would 
seem natural to evolve the position of a particle in terms 
of its X and y coordinates, without direct consideration of 
the third dimension (excepting the implicit information 
found in h{x, y)). 

This paper presents a 2D Langevin equation for diffu- 
sive particle motion over a curved Monge-gauge surface 
that may be static or fluctuating in time. The numeri- 
cal implementation of this equation leads to simulations 
of increased efficiency relative to the more general 3D 
random- walk methodology of Holyst et al. [TH] . Simula- 
tions are carried out here for surfaces with lateral peri- 
odicity. 

In the quenched limit, simulation results for a variety of 
differently-shaped surfaces (periodic in a box of linear di- 
mension L) suggest that the ratio between the projected 
and the curvilinear diffusion coefficients is equal to the 
ratio of projected membrane area, A± = L^, to actual 
curvilinear membrane area, A (to within the precision we 
are able to achieve numerically). That is, D/Dq ~ A±/A 
for static Monge-gauge surfaces. These results directly 
refute an earlier claim of non-integer scaling with surface 
area and disprove the suggestion [16 that the effec- 
tive medium approximation should be the most accurate 
approximation available for predicting D from Dq. 

To examine the effects of dynamic surface fluctuations, 
we consider Helfrich elastic membranes with thermal dy- 
namics described by a generalized continuum Langevin 
equation [S]. As a general trend, we show that the 
ratio D/Dq grows with increasing membrane fluctua- 
tion rates. A particle moving on a slowly fluctuating 
membrane exhibits D approaching the quenched disorder 
limit discussed above, while fast membrane fluctuations 
effect D values limiting to the somewhat larger annealed 
limit results discussed previously [16j |24]. For a mem- 
brane of given rigidity, the crossover between these two 
regimes is controlled by a dimensionless parameter x = 
kBT/{4riLDo), which represents the ratio between char- 
acteristic particle diffusion time and membrane shape- 
relaxation time scale in a hydrodynamic medium of vis- 
cosity Tj. A formal operator expansion in powers of 
demonstrates the pre-averaging Fokker-Planck equation 
[24] as the limiting behavior associated with fast mem- 
brane fluctuations (x~^ — > 0). Additionally, this expan- 
sion suggests higher-order corrections, which may be of 
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FIG. 1: A typical equilibrium configuration of a semiflexi- 
ble elastic membrane described in Monge gauge. Projected 
self-difltusion in the x ~ y plane (with effective diffusion co- 
efficient D) is necessarily reduced relative to the curvilinear 
self-diffusion over the actual surface (with diffusion coefficient 
Do) since any path followed on the membrane surface projects 
to a path of equal or smaller length in the base plane. 



use in further theoretical study. 

We consider a wide range of x values in our simulations 
to fully characterize the transition between quenched and 
annealed regimes. Typical membrane proteins in suffi- 
ciently large tensionless biomembranes are expected to 
exhibit intermediate x values, with D values falling in- 
between the quenched and annealed limits. 

II. LANGEVIN EQUATION FOR DIFFUSION 
ON CURVED SURFACES 

A. Membrane geometry 

Let us briefly review some elementary geometric con- 
cepts of surfaces in order to introduce notation [33] . A 
curved 2D surface is parameterized by two independent 
coordinates. In the Monge gauge, one makes use of two 
Cartesian coordinates {x, y) with respect to a flxed frame 
in a 3D Euclidean space; in the absence of overhangs [32], 
z = h(x, y) uniquely represents the membrane height 
from the x — y plane (see Fig. [l|. It is also convenient 
to deflne 0(r) = z — h{x, y), to yield the surface equation 
in 3D: 0(r) = where r = {x,y,z). The gradient vec- 
tor, Vr0 = {—dxh, —dyh, 1), points normal to the surface 
with magnitude equal to the determinant of the surface 
metric 

g^\Wr^\^ = l + {dxhf + {dyhf, (1) 

which formally determines the area rescaling of an in- 
finitesimal surface element upon projection onto the x — y 
plane, i.e. 

dA = ^dA± = Vgdxdy, (2) 
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where A represents the curvilinear area of the surface 
and A±, the area projected on the x — y plane. 

Using a condensed notation Xi ^ x,y for i — 1,2, the 
metric tensor elements of the surface in the Monge gauge 
read 



(3) 



with 6ij being the Kronecker delta and di = d/dxi [M] , 
Hence one obtains the determinant g = det(g) and the 
inverse tensor 



(g ^)tj^ % - 



dih djh 
9 



(4) 



Another useful quantity is the curvature tensor defined 
as [33] 



_ didjh _ d,dj(j) 
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(5) 



We will express some of our results in terms of the mean 
curvature, 



(6) 



which gives the average of the two inverse principal radii 
of curvature at a given point on the surface. The Einstein 
summation convention is assumed above and in all that 
follows. 



B. Discrete random walk on a membrane 

Our goal is to estabhsh a Langevin equation that de- 
scribes the Brownian motion of a particle confined to 
z = h{x,y). To this end, we first consider the case of 
discrete (but small) displacements over the membrane 
surface. Let us assume that a particle starts at Tq on 
the membrane surface (i.e., 0(ro) = 0) before the dis- 
crete jump. A Brownian jump of finite duration At 
may be broken down into two consecutive steps |18l 
(see Fig. [2]): i) The particle makes a random jump, 
— ({A(i},A^2), to a new position ri = tq -I- A^; 
this displacement is restricted to lie in the local tangent 
•plane to the surface at Tq and is distributed according to 
a Gaussian probability distribution to be specified below, 
ii) The final location, r(At), of the particle is then ob- 
tained by an instantaneous normal projection of ri back 
onto the membrane surface. 

The net displacement of the particle, Ar = r(Ai) — Tq, 
may thus be written as 



Ar = AC + A„ 



(7) 



where the second term on the right hand side represents 
the normal projection (step ii), i.e. 



v{M) - n 



(ri)V,,,^(ri) 
|Vr,0(ri)P 



(8) 
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FIG. 2: Schematic representation of Eq. ([7|. Every dis- 
crete random-walk step (of duration At) over a curved surface 
(shown here by a thick curve for simphcity) may be broken 
down into two consecutive steps [18) . First, a random jump of 
size and duration At from the initial position ro to a new 
position ri in the local tangent plane (shown by dashed line). 
Second, an instantaneous normal-projection move of size A„ 
from ri to the final location, r(At), on the surface. In the 
limit At 0, particle motion becomes exactly confined to 
the surface and, as shown in the text, may be describe by a 
well-defined Langevin equation. For the sake of representa- 
tion, the jump size in the figure is exaggerated; note that A„ 
is perpendicular to both dashed tangent lines within the first 
order in A<^. 



and for the first term we have 



AC-Vro0(ro)=O. 



(9) 



The above equations have been used to numerically sim- 
ulate diffusion over a curved surface [TH|. A significant 
merit of this methodology is that it may be used on 
arbitrarily shaped surfaces; i.e., a single-valued Monge- 
gauge parameterization is not required to implement this 
random-walk algorithm. In the case of a surface that may 
be specified by h{x,y), this procedure may be simpli- 
fied considerably, reducing the inherently 3D algorithm 
to a strictly 2D equation of motion, with the surface con- 
straint z = h(x,y) met exactly. (In contrast, the surface 
constraint is satisfied only to second order in jump size 
within the random-walk algorithm as implied by Eq. ([s]) 
[TH].) This 2D equation of motion is specified in the fol- 
lowing section. 

Let us consider the second term in Eq. ^. For suf- 
ficiently small At, the length of the initial jump, A^, 
becomes small enough to justify expanding 0(ri) around 
Tq. To second order in A^, 



m 

which shows that the projection displacement, A„, is of 
the second order in the random jump size, A^. 

Using the notation introduced in Section |II A[ and 
rewriting this equation in Monge gauge, one finds the 
following set of equations 



Ax,^-^(^AQAad,dkh^^ 



-AQ+0{Ae) (11) 
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Az = UAQACkdjdkh^ ^ 



AC. + 0{AC'). (12) 



Recall that k = 1,2 denote only the projected x — y 
components. 

Note that in the first term on the right hand side of 
both equations, only the x — y components of the random 
jump enter. By virtue of Eq. it follows that 

the z-component of the random jump, A^^ = AQdih, is 
not an independent variable. As a result, for sufficiently 
small time steps or random jumps, the displacement in 
the z direction, Az, is fully described by the projected 
displacement in x — y plane. In other words, Eq. ( |12[ ) 
does not constitute an independent equation but merely a 
relation for the particle height variation, Az, after each 
discrete step as enforced by the surface constraint z — 
h{x,y). Indeed, one can easily check using Eqs. (11) 
and ([9|, that Az is simply equal to the variation of the 
membrane height. Ah, for a small projected jump Axi, 



Az = Ah = Ax, a/i 



1 



Axi Ai 



■jd,djh\^ 



(13) 



where the right hand side is a Taylor expansion. Hence, 
in the limit of infinitesimal jumps, the particle is auto- 
matically confined to the surface as its x — y position 
evolves according to Eq. (111. 




C. Continuous time limit: Langevin equation 



FIG. 3: EfTect of curvature on the projected diflusive process: 
Curvature of a surface (here a line for simplicity) leads to de- 
terministic drift within the Langevin equation. Two random 
steps of equal magnitude (cl and cr) cover the same curvilin- 
ear distance irrespective of direction, but translate into dif- 
ferent projected distances {\xl — xo\ and {xn — xo\) due to 
surface shape. At a given position, the effect always favors 
displacement in the same direction, regardless of the exact 
(small) step size. Note that both slope and second derivative 
of the surface must be non-zero to break symmetry and drive 
the particle in a specific direction. 



We now consider the time evolution relation ( 1 1 ) in the 
limit of infinitesimal time steps, that is formally At —^ dt 
or for the random jumps, Ai^^ — > d^^. The random jump 
was chosen so that it locally describes a free Brownian 
motion in the local tangent plane (characterized by a 
Gaussian distribution function) with a curvilinear diffu- 
sion coefficient Dq and duration At. Thus, the projected 
jump probability density may be written as 



J dr Go(r) S{z - xAh) n,;<5(AC. - x^) 



: N exp 



/ dr Go(r) 5{z — Xidih) 
-5,jAC.AC,/(4Z?oAi) 



(14) 



where Go{x,y,z) = {'i-KDoAt)-^'^ ey.^[-v^ / AD^At] is 
the usual Green's function for homogeneous 3D diffusion 
and the normalization is given \yy M = ^/{AirDoAt). 
Formally, this distribution defines a generalized Wiener 
process ^35j . Ci(^)i (with the non-constant kernel gij) in 
the base plane, whose properties are well-defined in the 
infinitesimal time-step limit; in particular, one has 



dQdQ^2Do{g-% dt, 



(15) 



and d^i dt — 0. 

Rewriting Eq. ( 11 ) for an infinitesimal random jump, 
one thus obtains 



dxi ^ -Do 



(g ^)jkdjdkh 



d,h 



dt + dC» 



(16) 



Note that the second-order term in random jump (coming 
from the projection displacement A„) leads to a drift, 
which is of the order ~ dt, while the noise term scales as 
dCi ~ dt^/^ and all higher-order terms formally vanish. 

The corresponding Langevin equation immediately fol- 
lows as 



Xr{t) = 

where the drift term, Vi 



•H + Tij r]j{t), 
Vi[h{{xk{t)})], reads 



Vi = -Do 



(g ^)jkdjdkh 



dih „^ dih 
— = -2DoH —. 

9 ^/g 



(17) 



(18) 



It is induced by the mean curvature, H , and local gradi- 
ent of the surface. The fact that surface curvature drives 
local drift is easily understood in the context of diffusion 
on a ID curved line (see Fig. |3]) . In the presence of curva- 
ture and non-vanishing slope, equidistant displacements 
in opposite directions along the curve necessarily lead to 
projected displacements of different magnitude. 

In Eq. (17 1, we have defined Tijrjj^t) = dQ/dt, where 
rji (t) is a Gaussian-distributed white noise of zero mean, 
{rii{t)) = 0, and variance 

(77,(i) 77,(0) =2Do<5,j<5(t-0, (19) 

and Tij = Tij[h{{xk{t)})] is the square-root of the 
(positive-definite) inverse metric tensor, i.e. 

^ nkTjk- (20) 
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The noise term in Eq. (17) is thus multiplied by a factor 
containing explicit dependence on (x, y), a clear manifes- 
tation of "multiplicative noise" [U [3S] . The derivation we 
have presented suggests that the noise must be treated 
within the Ito prescription [JJ |3S]. The noise prefactor 
may be calculated by diagonalization of the inverse met- 
ric tensor as [551 



dihdjh 



(21) 



This completes the derivation of the Langevin equation 
for free diffusion on a curved surface in Monge gauge. 

The projected motion of a particle confined to a 
Monge-gauge surface may be regarded as formally equiv- 
alent to the motion of a particle in a 2D complex in- 
homogeneous environment characterized by an effective 
"potential", U{{xi}), and a locally varying diffusion ma- 
trix, 'Dij{{xk\), both determined by the surface metric. 
The latter may be read off directly from the noise term 
as 



(22) 



The explicit form of the effective potential may be seen 
more clearly from an equivalent Fokker-Planck equation 
[IS]. Using methods described in [H |33 [57], one can 
derive the following Fokker-Planck equation from the 



Langevin equation (17 1, i.e 



dtV{{x,},t) = Dod^ 



3 1 -^^J 





-11 
s/s). 


(23) 






(24) 


^)] 
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(25) 



which governs the time evolution of the projected prob- 
ability density, V{{xi\^t), with the normalization con- 
dition ^ d^xV{{xi\,t) = 1. In the last equation, Uj = 
—djU stands for the gradient of the effective dimension- 
less potential [IB] 



C/({xJ) = -^ln5. 



(26) 



III. DIFFUSION OVER A STATIC SURFACE 



Brownian Dynamics simulations on a quenched 
surface 



In order to demonstrate how the projected diffusion co- 
efficient, D, is affected by static surface undulations, we 
consider diffusion over periodic surfaces h{{xi + niL}) = 
h{{xi}) with i = 1, . . . , d in d = 1 and 2 surface dimen- 
sions and rii being integers (hence, in 2D, the surface con- 
sists of identical square-like patches of size LxL). This is 



a convenient choice for numerical simulations, and allows 
for comparison with existing works. 

The projected diffusion coefficient, D, is calculated by 
numerically evolving the projected position, Xi(t), of a 
single particle (or equivalently, many non-interacting par- 
ticles) according to the Langevin equation ( 17 1 and then 
evaluating 



D 



lim 

t-^oo 2dt 



Y,[^,{t)-xMV 



(27) 



where the bar sign denotes average over an ensemble of 
particle trajectories with random initial conditions, Xi{0) 
(as well as an ensemble of surface configurations in the 
case of the random surfaces introduced in Section III C 3[ ) . 
In all cases considered here, we choose surfaces (or ensem- 
bles of surfaces) with sufficient symmetry to insure that 
the principal directions x and y are equivalent. (The sum 
over i in Eq. ( |27[ ) only serves to enhance sampling.) 

In the simulations, we apply a standard iterative 
(Brownian Dynamics) algorithm [55] by discretizing the 
Langevin equation ( 17) using sufficiently small time steps 
of length St. Hence at time t = nSt (with integer n > 0), 



Xi{n+1) = Xi{n)+Vj[xk{n)]St+ ^/2D^Tij[xk{n)]w-i{n) , 

(28) 

where Wi are Gaussian-distributed random numbers 
with zero mean, {wi{n)) — 0, and unit variance 
{wi(n) Wj{m)) = 5ij Snm- This form of the discrete equa- 



tion relies on the Ito interpretation of Eq. ( 17 1 previously 
mentioned. We emphasize that the lateral position of the 
particle is not constrained to any discrete lattice. 

Within our level of description, the projected diffusion 
coefficient, D, must be linearly proportional to Dq. This 
is because the only possible dimensional quantities affect- 
ing D are Dq, L, and the surface shape, h{x,y), and the 
only possible dimensionless quantity incorporating Dq is 
the ratio D/Dq, which must be expressible in terms of 
the other dimensionless quantities in the physical prob- 
lem [Sg. That is, D/Do = f[h{x,y)/L] is a functional of 
the reduced surface shape alone (see also Section IV C 1 1 . 



In the case of surfaces described by a single parameter, 
A (such as the amplitude in the sinusoidal and cycloid 
surfaces specified by Eqs. (35) and (36) below), D/Dq is 
a function only of A/L. 

In the quenched simulations, we typically employed a 
rescaled time step of DoSt/L'^ ~ 10~^ — 10~^ and ran 
the simulations for approximately 10^ — 10^ steps. It 
was verified that this choice of time step leads to con- 
vergence of the reported results for the roughest surfaces 
considered; i.e., further reduction of time step did not 
change the results within reported error-bars. Statisti- 
cal averages are made over an ensemble of 5 x lO'^ — 10^ 
trajectories (and membrane configurations where neces- 
sary) as discussed above. The error-bars are calculated 
using standard block-averaging methods [40]. Note that 
for a typical protein in a lipid bilayer, the diffusion co- 
efficient Dq ^ 1 /zm^/s. Assuming a membrane patch 
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of size L ^ 0.1/zm, the simulation time steps fall in the 
range 5t ~ 0.01 — I^jls. 

We remark that simulations may also be done with the 
random-walk algorithm defined via Eqs. as dis- 

cussed in Ref. [18]. As we have shown in the previous 
Section, for sufficiently small time steps, this method is 
equivalent to the Langevin approach. The latter has the 
advantage that it implicitly constrains the particles to the 
membrane. Furthermore, it is a dynamical equation that 
may be generalized to study more complicated systems 
with, for instance, inter-particle or particle-membrane 
interactions. The random-walk algorithm is more con- 
venient for surfaces with complicated geometry (such as 
3D supramolecular assemblies) ,18, ,41] , where the present 
Monge-gauge Langevin equation is not adequate. Also 
for surfaces with very large gradients, it is clear that the 
Langevin approach will eventually require smaller time 
steps than the 3D random-walk approach and thus longer 
computational runs. In practice, however, for the exam- 
ples considered below, the Langevin simulations appear 
to be faster (typically by about 50% computer time in 
two dimensions) as comapred with the random-walk al- 
gorithm. There are simply less computations per time 
step in the Langevin approach due to the analytical re- 
duction of complexity when a Monge-gauge description 
is possible. Random-walk simulations were carried out 
on a number of the reported test-cases for comparison 
purposes-always with identical results to the Langevin 
simulations. For simplicity, these random- walk results 
are explicitly displayed only in Figure [4] 

We will first consider numerics on simple periodic sur- 
faces in one dimension where analytical results are avail- 
able. We then progress to simple 2D surfaces (i.e., a 
sinusoidal surface and a crested cycloid j23j) and finally 
the case of an elastic membrane. 



B. One dimension: exact result and simulations 



In ID, the local diffusion matrix (22 1 is a scalar and 
reads V = Do/g. Highly steep surface segments with 
large slopes relative to the x-axis {g ^ 1) represent re- 
gions with small local diffusion coefficient, 'D/Dq ^ 1, 
and large negative effective (dimensionless) potential, 
—U^l, Eq. (26). Regions of high slope may be re- 
garded as trapping potential wells (see also Fig. [sjbelow). 
Physically, the cause of this effect is clear. In regions of 
high slope, a particle which travels a large distance along 
the contour of the curve travels a relatively much smaller 
projected distance. The majority of travel is occurring 
in the z direction, not contributing to the projected dif- 
fusion coefficient, D. 

The exact ID expression for the projected diffusion co- 
efficient along an arbitrary periodic curve may be written 
in terms of the metric determinant, g = 1 + (dxh)^, as 




FIG. 4: Rescaled projected diffusion coefficient as a function 
of rescafed curve amplitude for free particfe diffusion over a 
quenched ID curve defined by Eq. ( 30 1 for two values of B = 
and B = A as indicated on the graph. Symbofs show sim- 
ufation data using Langevin dynamics (open diamonds and 
circles) and random-walk algorithm [18] (cross an d p lus sym- 
bols). The solid curves show the exact result, Eq. (29 1. Error- 
bars are smaller than the size of symbols. 
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where brackets denote a contour average defined as 
(...) = J^dx^(...)/ J^dxy^, and Lc is the contour 
length of the curve over a single projected period L. 

It is not surprising that analytical results are available 
in ID. The usual ID Green's function for diffusion pre- 
dicts statistically Gaussian behavior for motion along the 
contour length of the curve. For a Monge-gauge curve, 
this solution is readily transformed to the x-axis pro- 
jected results, since contour length and projected length 
may be unambiguously mapped to one another. The pro- 
jected solution in turn leads to Gaussian behavior in the 
long time limit with a simple rescaling by the square of 
the projected length to contour length ratio [17] . 

In Figure [4] we plot the Langevin simulations results 
(open symbols) for a quenched periodic curve defined as 

, / s . 27ra; „ Anx 

h{x) = A cos — — -I- B cos — — . (30) 

Lj Li 

for two different values of S = and B — A. As seen, 
the numerical results exactly coincide with the analyti- 
cal curve, Eq. (29 1, confirming thus the validity of the 
present Langevin approach. For comparison, we have 
also plotted in the Figure the simulation data (cross and 
plus symbols) obtained by applying the random-walk al- 
gorithm [18J, Eqs. ([7])-([9|, which, as expected, coincide 
with both Langevin simulations data and exact results. 



12 3 4 

AIL 



FIG. 5: Rescaled projected diffusion coefficient as a function 
of rescaled surface amplitude for free particle diffusion over a 
quenched surface defined by Eq. ( 35 1 for two values of B = 
and B = ^ as indicated on the graph. Symbols show Langevin 
simulation data (open diamonds and circles), while solid and 
dotted curves show the area-scaling and the effective medium 
predictions, Eqs. (34 1 and (331, respectively. Inset shows 
A=l and L = 2-k. 



surface profile for 5 = 0, 
smaller than or equal to symbol size. 



Error-bars are 



C. Diffusion on a 2D quenched manifold 

Unlike the ID problem, a rigorous result for the pro- 
jected diffusion coefficient is not known for a quenched 
2D surface. This relates to an intrinsic difference between 
Brownian motion on ID and 2D manifolds ^16]. Namely, 
on a ID manifold, the distance between two points is 
given by the contour length which provides a natural 
parametrization, while on a 2D manifold, many differ- 
ent paths connect two given points, leading to analyti- 
cal difficulties for establishing an exact result for a static 
surface [HIIT^. Nonetheless, several approximate results 
have been presented in the past that we shall briefly re- 
view below following the discussions by Gustafsson et al. 



1. Approximate analytical results 

First, one may establish an upper bound and a lower 
bound that limit the admissible values for D using vari- 
ational considerations (THJ HI] . The result for the upper 
bound in the quenched surface limit reads |16j 



D_ 
Do 



1 



upper 



(31) 



where brackets denote a surface average defined as 

(...) = /„^d2xV5(...)//„VxV5. 

Intuitively, one expects that the quenched ID result, 
Eq. ( p9| , establishes a lower bound since correlation ef- 
fects induced by quenched surface undulations are ex- 



FIG. 6: Same as Fig. [5] but for B — 0. Symbols represent 
simulation data, dashed and dot-dashed curves are the pre- 



dicted upper and lower bounds, Eqs. (31 1 and (32 I, and solid 



and dotted curves show the area-scaling prediction, Eq. ( 34 1 



and the effective medium prediction, Eq. (331, respectively. 



pected to be stronger in ID [T7] . Using variational meth- 
ods, an expression similar to the ID result may be ob- 
tained as [T3] 



D 



lower 



(l/.9> 



1 



(32) 



In addition to these variational bounds, a mean-field 
approach-referred to as effective medium approximation 
(EMA)-has been introduced [121 EI], which is supposed 
to provide a good approximation for projected diffusion 
coefficients on strongly disordered quenched surfaces, i.e. 



D 



1 



EMA 



(33) 



Finally, another plausible approximation for D can be 
deduced based on naive scaling arguments. On dimen- 
sional grounds, one might expect that an effective local 
isotropic projected diffusion coefficient would scale with 
the ratio of local projected area to actual curvilinear area, 
i.e., D/Do ~ d^_L/d^ = (see Eq. Averaging 

this result over the entire surface leads to 



D 



area scaling 



1 

7^ 



A 



(34) 



where the second equality follows from the definition of 
the averaging process as used in Eq. (31). One could 



have also guessed the second equality by applying similar 
arguments to the entire surface area instead of making 
local considerations. This conjecture finds some support 
in the fact that it is exact in ID (see Eq. (29 1). To our 
knowledge, no systematic derivation has been given so 
far in favor of the area-scaling prediction, though it has 
been discussed in several previous works [TSJH3J|2S|. 
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FIG. 7: Rescaled projected diffusion coefficient as a function 
of rescaled surface amplitude for free particle diffusion over a 
quenched crested cycloid defined by Eq. (36 1. Symbols show 
Langevin simulation data (open squares), and the solid and 
dotted curves show the area-scaling and the effective medium 
predictions, Eqs. ( |34[ | and ( |33[ ), respectively. Inset shows the 
height profile of a cycloid with L — 2tt and A = 0.7. Error- 
bars are smaller than or equal to symbol size. 



2. Simple periodic surfaces 



In order to examine the aforementioned analytical ap- 
proximations, we first consider a simple periodic surface 
defined as 



in a parametric form as 
A 

A 



2ttu 



, sm 

2tt J \ L J 



A^ f2TTU 

cos 

2ttL V L 



( 2t:v\ 



(36) 



where z = h{x{u, v), y{u, v)) represents the surface height 
(see Fig. [t] inset) and u and v are real numbers. Note 
that the cycloid is not single- valued for A/L > 1 and at 
A/L = 1, it shows singular cusp-like peaks. Therefore, 
we restrict our discussion to the range with A/L < 1. 
Diffusion on this surface was previously investigated in 
Ref. [23], where it was suggested that the projected diffu- 
sion coefficient deviates from the area-scaling prediction 
and follows the scaling law D/Dq = {A±/Ay''^'^. 

In Figure |7] we have plotted the results from our simu- 
lations and compared with both the area-scaling and the 
effective medium approximation (solid and dotted curves, 
respectively) . It is evident that the data closely follow the 
simple area-scaling rule D/Dq = {A±/A). In the present 
case, the area-scaling and effective medium curves di- 
verge for growing A/L ratio, with the latter systemati- 
cally underestimating the projected diffusion coefficient. 
The spurious results presented in Ref. [23] apparently 
stem from an incorrect treatment (neglect) of the surface 
metric within their simulation protocol. 



Quenched elastic membrane 



2ttx 2iTy 

hix, v) — Acqs cos 

^ L L 



Anx Any 

B cos cos . (35) 

L L 



Using Brownian Dynamics simulations as described in 
Section [III A| we calculate the projected diffusion coef- 
ficient for two different cases with B = and B ~ A 
(shown by open symbols in Fig. [5|. As seen, in both 
cases the simulation data for D/Dq as a function of A/L 
agree well with the area-scaling prediction (solid curves) . 
There are strong deviations from both upper and lower 
bounds as shown in Figure |6] (compare with dashed and 
dot-dashed curves, respectively). But for the above sur- 
face, the area-scaling result stays very close to the ef- 
fective medium prediction (dotted curve). Note that for 
small A/L, the surfaces are nearly flat and different ap- 
proximations work similarly well in this limit. Only in 
the large A/L limit, can one distinguish between differ- 
ent predictions. For A/L ^ 1, the upper bound value 
for D/Dq tends to 1/2, while the lower bound decays as 
^{A/L)~'^. The simulation data (as well as area-scaling 
and EMA predictions) fall off as --(A/L)~^ 

Next we consider diffusion on a crested cycloid defined 



We now turn our attention to the case of a randomly 
undulated, quenched elastic membrane of size L x L 
with periodic boundary conditions. In the continuum 
limit (large length scales), we use the Helfrich membrane 
model with the Hamiltonian 



n=^ [ dp [K{vlh)' + aiyphf] 

^ J A 



(37) 



2L2 



where p = (x,y), K is the bending rigidity, a the surface 
tension, and 



(38) 



the energy spectrum of membrane fiuctuations. For sim- 
plicity, we shall focus here only on two special cases, 
where only one of the elastic coefficients is non-zero. 

At equilibrium, any configuration of the membrane 
may be viewed as a linear combination of several Fourier 
modes, q — {qx,qy), similar to Eq. (351 but with ampli- 



tudes, ft.q, distributed according to the Boltzmann weight 
~ exp(— Ti/fcB?^)- For numerical purposes, one needs 



9 




1 


"1 1 

|7-' /■' 


1 1 ' 

y 


o 


• ■ / ' 

["^^•'' /' 








0.1 


/ 0.2 

/ -0-2 


"^•> ■ 




/ 1 


^ . ♦ ---^ y 






0.01 0.1 


1 10 



FIG. 8: Typical height profile, h{x,y), metric determinant 
profile, g{x,y), and effective (dimensionless) potential pro- 
file, U{x,y), for an elastic membrane with bending rigidity 
K/{KbT) — 0.01 and surface tension a — along the line 
y = 0. Vertical lines show typical regions with minimum (thin 
lines) and maximum (thick lines) potentials. As discussed in 
the text, these curves reflect a simple ID slice through the 
surface. Consequently, non-perfect correlations are observed 
between, e.g., extrema in h{x,0), g{x,0) and U{x,0), reflect- 
ing the possibility of j/-component contribution. 



FIG. 9: Rescaled projected diffusion coefflcient, D/Do, of a 
freely diffusing particle on a quenched ruffled membrane as a 
function of the rescaled membrane bending rigidity K/ (fceT) 
for M = 32 and surface tension a — 0. Symbols represent 
simulation data, the dashed curve is the upper bound, Eq. 
(311, the solid curve is the area-scaling prediction, Eq. 
the dotted curve is the effective medium result, Eq. 



and the dot-dashed curve shows the predicted lower bound, 
Eq. ( |32[ ). Inset shows a typical membrane configuration with 
K/{kBT) = 0.01, L = 1 and M = 32. Error-bars, ff not 
shown, are smaller than or equal to the size of symbols. 



to introduce a small-scale cut-off, a (mimicking, for in- 
stance, a specific microscopic length scale) in order to 
filter the irrelevant short-wave- length modes. The height 
profile is still continuous and defined in space via 



Hp) 



(39) 



where q = {qx,qy) = (2Tm,2T:m)/ L with n,m being in- 
tegral numbers in the range —M/2 < n,m < M/2 and 
M ~ L/a. In order to investigate diffusion over elastic 
membrane surfaces, we first numerically generate many 
equilibrium membrane configurations consistent with the 
above Hamiltonian by drawing hq from normal distribu- 



tions for each mode [121 IB] (see also Section IV I 

We simulate the diffusion of a particle using the 
method described in Section III A and average the mean- 
square displacement over an ensemble of 5 x 10'^ — 10'' 
different equilibrium membrane configurations. As pre- 
viously discussed, the rescaled projected diffusion coef- 
ficient, D/Dq, is only a function of dimensionless mem- 
brane shape parameters, which, in the present case, are 
the rescaled bending rigidity K/ik-^T), rescaled surface 
tension cra^ / {k^T) and the ratio M — Lja (see also Sec- 
tion IV C 1 1 . Here we choose M = 32 for numerical 
convenience. (This choice is most sensible assuming a 
membrane patch of size L ~ 100 nm, which implies a 
small-scale cut-off of a ~ 3 nm, comparable to the size of 
lipid molecules.) 

A typical membrane height profile along the y = axis 
is shown for cr = and K/{k-QT) = 0.01 in Figure [s] to- 
gether with the profile of the surface metric determinant 
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FIG. 10: Rescaled projected diffusion coefficient, D/Dq, of a 
freely diffusing particle on a quenched ruffled membrane as a 
function of the rescaled membrane surface tension aaF' / {ksT) 
for M — 32 and bending rigidity K = 0. Symbols represent 
simulation data, the dashed curve is the upper bound, Eq. 
(311, the solid curve is the area-scaling prediction, Eq. (34l, 
the dotted curve is the effective medium result, Eq. (33 1, 



and the dot-dashed curve shows the predicted lower bound, 
Eq. (32 1. Inset shows a typical membrane configuration with 
aa^/(fcSr) = 0.2, L = 1 and M = 32. Error-bars, if not 
shown, are smaller than or equal to the size of symbols. 



and the effective dimensionless potential, Eq. (26 1, ex- 
perienced by a diffusing particle. Potential wells (peaks) 
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correspond to steep (flat) segments of the membrane and 
therefore a large (small) metric determinant as indicated 
by vertical lines on the graph. (Note that both x and y 
components of the surface gradient contribute to g{x, 0) 
and t/(a;,0), and that only the x-component can explic- 
itly be seen in h{x, 0) plot in the Figure.) 

The results for the disorder- averaged projected diffu- 
sion coefficient are plotted in Figure[9]as a function of the 
rescaled bending rigidity at fixed ct = (open symbols) 
and in Figure [lO] as a function of the rescaled surface 
tension at fixed K — Q. Both figures show qualitatively 
similar behavior for the simulation data and theoretical 
predictions. (It is to be noted that the theoretical curves 
are computed using Eqs. (31 )-(34l and directly from the 



same ensemble of membrane configurations that are used 
in our Langevin simulations.) For large bending rigidi- 
ties (or large surface tensions), the membrane stays close 
to a flat configuration and the data as well as all other 
theoretical predictions for D/Dq tend to one (free planar 
diffusion) in the same manner. The difference between 
these results becomes pronounced and they can be dis- 
tinguished in the limit of small bending rigidity (or small 
surface tension). 

Interestingly, we find that in both cases the simula- 
tion data agree well with the area-scaling prediction ( 34 ) 



(shown by a solid curve) down to very small values of 
K and cr. In agreement with previous examples studied 
in this paper, the effective medium approximation (dot- 
ted curve) diverges from the simulation data for small K 
(or small a) and, despite the suggestions made in Ref. 
|16j . it represents a poorer approximation as compared 
to simple area scaling. 



IV. DIFFUSION OVER A DYNAMICALLY 
FLUCTUATING SURFACE 



namic medium is customarily described by a Langevin- 
type equation [8tl44ll45]. 



dh{p,t) 
dt 



/ + 00 
dp'A{p^p')[F{p\t) + r^{p\t)], (40) 
-oo 



where A(p— p') is the hydrodynamic interaction specified 
by the Oseen tensor (a no-slip boundary condition on the 
impermeable membrane surface is assumed) j45l |46l |47] , 



Hp-p') = 



8TTr]\p- p'\ 



(41) 



F{p, t) is the force per unit area given by the functional 
derivative of the Hamiltonian, 



F{p,t) 



sn 



Sh{p,t) 



(42) 



and ri{p, t) is a spatially-varying Gaussian, white ran- 
dom force of zero mean, {'r]{p, t)) — 0, that satisfies the 
fluctuation-dissipation theorem 

(77(p, t) rjip', t')) = 2 k^TK-\p - p') 5{t - t'). (43) 

This choice of noise clearly ensures that the membrane 
relaxes to equilibrium at long times. 

To proceed, we express the membrane Langevin equa- 



tion (40 1 in terms of Fourier modes /iq(i) with a discrete 
spectrum of wave- vectors q = (qx,qy) = (2nn,2Trm) / L 
where —M/2 < n,m < M/2 and M = L/a as discussed 
in Section pIIC3[ Hence, we have 



(44) 



using which we may write Eq. (40 1 as 



So far we have only considered self-diffusion on a 
quenched surface. In this Section, we proceed by inves- 
tigating dynamical effects from temporal surface fluctu- 
ations on the particle diffusion over a two-dimensional 
Monge-gauge membrane. As it turns out, the effects of 
membrane fluctuations on the projected diffusion coeffi- 
cient can vary considerably depending on how fast mem- 
brane undulations relax relative to particle diffusion over 
the surface. In order to quantify this observation, one 
first needs to specify the dynamics of the surface. 



A. Membrane dynamics 

We shall adopt the Helfrich membrane model described 
by a continuum time-dependent height profile, h(p,t) 
with p — {x,y), whose equilibrium behavior may be de- 
termined by a Hamiltonian T-C[h{p)]. As before, we as- 
sume periodic boundary conditions over a length scale 
of L in lateral directions. The dynamics of membrane 
height fluctuations in a low-Reynolds-number hydrody- 



(45) 



where Fq{t) is the Fourier-transform of the force (42 1, 
Aq = 1/ (477 q) is that of the Oseen interaction and ^^t) 
is a redefined Gaussian white noise with (^q(t)) — and 



(eq(t)eq'(i')>=2^q,-q' S{t~t'). 



(46) 



It is understood that because h{p,t) is real, hq{t) is a 
complex number with /i* (t) — ft._q(i), guaranteeing that 
the number of independent Fourier modes remains . 
Equation (45 1 introduces an efficient algorithm, known 



as Fourier-Space Brownian Dynamics method [421 146] . 
for simulating dynamics of biological membranes with 
arbitrary interaction potentials. 



B. Coupled membrane-particle equations 

In what follows, we shall focus on the case of a free 
elastic membrane described by the elastic Hamiltonian 
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(37 1. The force term in Eq. ( |45[ ) thus simphfies to a 
hnear form Fq(t) = —flqhq(t), with ilq defined in Eq. 



(381 



Putting together Eqs. (17 1 and (45), we have the fol- 



lowing set of dynamical equations for free particle diffu- 
sion on a free Monge-gauge elastic membrane, 

ii{t) = + Tij rij{t), (47) 
hq{t) = -ojq hq{t) + ^kuTL^Aq eq(t) , (48) 

where the decay rates, Wq, are defined as 

LUq = Aq r!q = Aq {K + Uq^) (49) 

It is to be noted that the particle dynamics (via Eq. 



(47l) in the present description is passive in the sense 



that it does not infiuence dynamics of the membrane (see 
the Discussion). However, as the height profile modes 
evolve in time (via Eq. (48 1), they directly infiuence the 
curvature-induced drift, Vi, and the multiplicative noise 
factor, Tij, which are related to the local height of the 
membrane in real space, h{x, y, t), as discussed in Section 

El 



The free membrane equations ( 48 ) can be solved an- 



alytically as they correspond to a set of independent 
Ornstein-Uhlenbeck processes [?]. The time evolution of 
hq{t) (starting from the initial value = hq{0)), 
erned by the conditional probability distribution 



P{hq,t\hlO) 



X exp 



2L2 



1 - e- 



.(50) 



C. Crossover from quenched to annealed disorder 
limit 



1. Reduced representation 

In order to demonstrate this feature, we shall use 
a reduced (dimensionless) representation by rescaling 
length and time scales by a characteristic length and 
time, that is, using x = x/L, or equivalently, q = qL, 
t = Dot/L^, as well as h{x,y,t) — h{x,y,t) / L, or equiv- 
alently, /iq(t) = hq{t)/L^. Hence we obtain 



Xi(i) 



hqii) = -Xll'q/lq(t) 



_ (51) 
Xh^^t), (52) 



where the rcscaled functions, namely, rescaled drift, Vi, 
rescaled multiplicative noise factor, fij, as well as rescaled 
noise terms fjiii) and whose explicit forms are given 

in Appendix [X] are independent of system parameters. 
The rescaled decay rate reads iuq — Aq fiq where f2q — 
kq^ + aq^ and Aq = 1 /q. 

Note that all physical system parameters have been 
condensed into four dimensionless combinations. These 
are the rescaled bending rigidity K ~ K/{kBT), rescaled 
surface tension a — aL^ /{k^T), M — L/a, and an addi- 
tional parameter 



X 



kBT 
4ijDoL ' 



(53) 



which we shall refer to as the dynamical coupling parame- 
ter as it embeds in itself the dynamical coupling between 
equations (47 1 and (48 1. Indeed, for vanishing dynami- 
cal coupling parameter x 0, the membrane equation 
reduces to hq{t) = 0, representing a static surface. 

One may also examine the conditional probability in 
rescaled units 



p(hq,m,o) = 




t\2 



,(54) 



In general, the coupled dynamical equations (47) and 



( 48 1 can be used to study the diffusion process on a 



surface that may be static, as studied in previous Sec- 
tions, or dynamic. The quenched disorder as well as the 
annealed disorder limits can thus be realized as special 
cases of Eqs. (|47| and (48), by choosing the membrane- 



which is defined via P(/iq, , 0) = P{hq,t\h°q,0) . 
For X — ^ 0, the static surface solution is recovered, that is 
P{hq, ilh^, 0) S{hq — h^), while for x ^ oo, it rapidly 
decays to the equilibrium distribution 



related parameters (such as bending rigidity or medium Pihq,t\hq,0) Poq(/iq) 



viscosity) such that the membrane decay rates, Wq, tend 
to zero (static surface) or infinity (fast-fiuctuating sur- 
face). Physically, however, these same limits can also 
be achieved by taking a very large or small curvilin- 
ear particle diffusion coefficient, Dq. The coupling be- 
tween particle and membrane motions introduces a new 
dimensionless parameter to the problem not present in 
the quenched case (or the annealed case). This parame- 
ter can be used to track the transition from quenched to 
annealed limits. 



27r 



1/2 



-r!q|/lq|V2 



(55) 



These two limits therefore refiect the static surface 
(quenched disorder) and the fast-fiuctuating surface (an- 
nealed disorder) limits, respectively. Solving the set of 
equations (51 1 and (52 1 for a finite intermediate value of 



X allows one to investigate the intermediate regime, as 
we shall do later. 

The reduced representation explicitly demonstrates 
that any dimensionless quantity obtained via simula- 
tions must depend on individual system parameters only 
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through the dimensionless combinations K, a, Xi ^'^d 
M = L/a. This imphes that the ratio of the projected 
diffusion coefficient and the curvihnear diffusion coeffi- 
cient, D/D{), has the general form 



D 



(56) 



which could be surmised simply from dimensional con- 
siderations [39] , From a practical perspective, the re- 
duced approach simplifies our investigation of the de- 
pendence of D/Dq on various system parameters. The 
addition of membrane motion adds only a single dimen- 
sionless parameter beyond those necessary in describing 
the quenched limit. This parameter, may be thought 
of as a rescaled solvent viscoscity. Since x appears only 
when surface dynamics is important, we tend to focus 
attention on it when calculating D/Dq over fluctuating 
surfaces. However, we stress that all the dimensionless 
parameters of the system influence the value of D/Dq. 

Another advantage of the reduced representation is 
that it immediately suggests a "small parameter" for 
use in a purturbative expansion around the analytical 
annealed limit results. In Appendix [C] we present an 
operator expansion in powers of x~^- zeroth order, 
the pre-averaging Fokker-Planck equation of Reister and 
Seifert [25 is recovered. Higher-order corrections are pre- 
dicted, but not pursued in this work. 



2. Characteristic time scales 



bending rigidity exhibits quenched behavior at long wave- 
lengths. 

It is evident that in a finite or periodic system, such as 
the ones considered in this paper where the wave- vector is 
bounded from below by the surface periodicity q > 2'k/L, 
the transition from one limit to another can be estab- 
lished theoretically by changing the dynamical coupling 
parameter over the range < % < cxd at fixed bending 
rigidity, surface tension and M. This is indicated by the 
second equality in Eq. (57 1. In what follows, we shall 
focus only on the case with zero surface tension, a = 0, 
where quenched behavior can be achieved within reason- 
able values for system parameters. 



D. Coupled particle-membrane simulation method 

In order to simulate diffusion on a dynamic surface in 
an arbitrary dynamical situation, we apply the standard 
Brownian Dynamics algorithm by discretizing the time 
into steps of sufficiently small length similar to what we 
described in Section |III A[ In order to bring up the role 
of the dynamical coupling and provide a unified analysis 
of the quenched-annealed crossover behavior, it is most 
convenient to simulate the system in the reduced rep- 
resentation using Eqs. (51 1 and (52 1, where all various 



system parameters condense into a few relevant dimen- 
sionless parameters. 

We evolve particle position using small rescaled time 
steps St such that at time t = ndt (with n > being an 
integer) 



The physical meaning of the dynamical coupling pa- 
rameter, Xi may be appreciated further when looking at 
the typical time scales that characterize the dynamics of 
particle diffusion and that of membrane relaxation. 



Xi{n+ 1) = Xi{n) + Vi{n) &i + \ 25iTij{n) Wj{n) (58) 

for i,j — 1, 2, where Wi are Gaussian-distributed random 
numbers with zero mean, (^^(n)) = 0, and unit vari- 



The typical diffusion time over a length of size q ^ may ance {wi (n) Wj (m) ) 



Other rescaled functions 



be written as = l/{Doq ), while the relaxation time of 
an undulation of wave-length q~^ is given by = uj~^. 
Using the definition of the decay rate, Wq, we have the 
time scale ratio 



Kq^ 



irjDoq 



Kq^ 



(57) 



In general, one can define the quenched disorder limit as 
the limit where the membrane relaxation time becomes 
much larger than the diffusion time. While, the converse 
limit of vanishingly small relaxation time represents the 
annealed disorder limit. When the two time scales are 
comparable, one encounters an intermediate situation. 

Clearly, at short wave-lengths q oo, one always deals 
with an annealed situation. But it is the long wave- 
length, g — !■ 0, behavior which is of primary practical 
importance in the context of physical measurements. A 
membrane with a finite surface tension always exhibits 
annealed behavior at sufficiently long wave-lengths as 
seen from Eq. (57), whereas a membrane with only 



appearing in the above equation are given explicitly in 
Appendix |A] They are evaluated from the instantaneous 
membrane height ft.({xj(n)}, n) — h(^x{n),y{n),rij ob- 
tained from the Fourier modes hq{n), via Eq. (44 1. 



Note that to evaluate hq at each time step, one can 



either use the discretized form of Eq. (52 1, i.e 



hq{n + 1) = hq{n) - x^^qh^^n) 6t + ^ xAq (5t Rq{n), 

(59) 

(with random complex variables i?q of zero mean and 
variance {Rq{n) Rq' (m)) = 2(5q,_q' (5„„i)i oi' directly 
draw hqin) values from the Ornstein-Uhlenbeck Gaus- 
sian distribution, Eq. ( [54| ). The second method has the 
clear advantage that it is based on an exact result and is 
computationally easier to implement as previously per- 
formed in Ref. [42]. We have used both methods in our 
simulations with the same results obtained within the 
statistical error-bars. 

The time step 5t must be chosen such that it ensures 
the convergence of Eq. ([58| (and Eq. ( 59 1 , if applied) . 
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Even if the membrane dynamics are drawn from the exact 
distribution ( 54 1 , care must still be taken to insure that 



a small enougOi is used to properly capture the effect of 
membrane dynamics in Eq. ( [58| . Both the drift and dif- 
fusion tensor for the particle depend upon the local height 
of the membrane and it is to be expected that at large 
X, one needs to use time steps limited by the membrane 
dynamics as opposed to particle diffusion. We typically 
choose rescaled time steps in the range 5t = 10^^ — 10"^, 
and run the simulations for 10^ — 10^ time steps. (As 
indicated previously, step size is verified by checking for 
convergence of results.) In the simulations, we measure 
the projected diffusion coefficient via the mean square 
displacement of a Brownian particle as defined in Eq. 
(27 1. (Simulated mean square displacements converge 



more quickly with time to a diffusive long-time behavior 
for large x or highly fluctuating membranes.) We average 
the results over an ensemble of 5 x 10'^ — 10'' particle tra- 
jectories with random initial conditions for the particle 
as well as the membrane height profile. 

In our simulations, the coupling parameter is varied 
over the range x = 10^* — 10^ and the rescaled mem- 
brane bending rigidity in the range K = 0.01 — 10, with 
M = L/a = 32 being chosen as in the quenched study 
Recall that at fixed bending rigidity 
kBT/{4:?]DoL) may be 



in Section |IIIC3| 
and system size, changing x 



viewed as changing either the curvilinear diffusion coef- 
ficient, Dq, or the medium viscosity, 77, or a combination 
of both. They all have the same quantitative effect on 
D/Dq, once K and M are fixed. 
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FIG. 11: Rescaled projected diffusion coefficient, D/Do, of 
a freely diffusing particle on a dynamically fluctuating elas- 
tic membrane as a function of the rescaled membrane bending 
rigidity K/ (fceT) for Af = 32 and surface tension cr = 0. Sym- 
bols represent simulation data for various dynamical coupling 
parameters, Xi spanning the whole intermediate dynamical 
coupling regime from quenched membrane limit with x = 
(filled squared, same data as in Figure |9| to a rapidly fluc- 
tuating membrane with x = 10 (open triangle- ups) . The 
dashed curve is the annealed theoretical prediction (as calcu- 
lated in Appendix [b|, while the solid curve is the quenched 
area-scaling prediction, Eq. ( |34[ |. Dotted lines are guide to 
eyes. Error-bars are removed for clarity (they are typically 
about the symbol size and comparable to or smaller than those 
shown in Fig. [9|. 



E. Simulation results 

In Figure [11] we show simulation results for the ratio of 
the projected diffusion coefficient to the curvilinear dif- 
fusion coefficient, D/Dq, as a function of the membrane 
bending rigidity with a number of different values cho- 
sen for the dynamical coupling parameter, x- The case 
of X = (flUed squares) represents a strictly quenched 



dashed curve in Figure 11 1 exhibiting a reasonably good 



membrane, with results reproduced from Section IIIC3 



As mentioned before, the quenched data agree with the 
area-scaling prediction (solid curve). (In the Figure, we 
use a log-log plot and explore bending rigidities down to 
very small (unphysical for lipid bilayers) values to clearly 
display the predicted trends. At large K, both quenched 
and annealed results converge to within just a few per- 
cent.) 

It is clearly seen that, as the dynamical coupling pa- 
rameter increases, the projected diffusion coefficient de- 
viates from and becomes larger than its quenched value, 
indicating that dynamical fluctuations of the surface en- 
hance particle diffusion. This is to be expected based 
on very general considerations that imply fast surface 
fluctuations help to excite the particle through potential 
barriers more easily as it moves over the surface [16] . For 
larger values of x (stronger surface fluctuations) , the data 
tend to the theoretical annealed prediction (shown by a 



quantitative agreement for x > 10 (see below). The cal- 
culation of the annealed curve is presented in Appendix 
[B] It slightly differs from the approximate analytical re- 
sults previously obtained in Refs. [TSJ [H] as we explain 
in the Appendix. 

In order to examine the crossover behavior from the 
quenched to annealed disorder limit, we show in Figure 
[T2I the simulated diffusion coefficients as a function of 
the dynamical coupling parameter at a flxed value of the 
bending rigidity K/ik^T) = 0.1. Here we change x over 
a range of six orders of magnitude spanning the transition 
from quenched to annealed regimes. The data exhibit a 
smooth crossover from the quenched area-scaling value 
(shown by a solid horizontal line) to the predicted an- 
nealed value (shown by a dashed horizontal line). The 
agreement is almost perfect at both converging ends 
within the statistical error-bars of the data (compara- 
ble to the symbol size, amounting to relative errors of 
about 2%). Note that in these regions of x (^-g-, for 
X > 1 or X < 10^^), the data exhibit a very weak de- 
pendence on X- As previously mentioned, this directly 
translates to a weak dependence of the ratio D/Dq on 
the curvilinear diffusion coefficient, Dq, when all other 
parameters are flxed (see Eq. (56)). In other words, in 
either limit x — or x — > 00, D/Dq depends only upon 
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FIG. 12: Rescaled projected diffusion coefficient, D/Dq, of a 
freely diffusing particle on a dynamically fluctuating elastic 
membrane as a function of the dynamical coupling param- 
eter X = ^bT/ (ArjDoL) at fixed membrane bending rigidity 
of K/iksT) = 0.1 for M = 32 and surface tension a = 0. 
Symbols represent simulation data, the dashed horizontal line 
specifies the corresponding theoretical annealed value (Ap- 
pendix |b]|, while the solid line represents the quenched area- 
scaling value (dotted curve is guide to eyes). A smooth 
crossover is observed between the asymptotic annealed and 
quenched predictions. Error-bars are comparable to the sym- 
bol size. 



the dimensionless parameters, K/{k^T) and M (as well 
as aL^ /{k^T) when present). The practical definition of 



these limits is implied by the curve in Figure 12 



We find the same trend for all values of the bending 

Qualitatively, 



rigidity as may be seen from Figure 11 



one can think of the x value where D/Dq reaches the 
mean value of the quenched and annealed predictions as 
a conventionally chosen crossover point beyond or be- 
low which the system is dominated by either quenched 
or annealed features as described in Section IIV C 21 In 
Figure [12] this corresponds to a value of x ~ O-l- For 
smaller bending rigidities, this crossover point tends to 
larger values of x (as may already be seen from Figure 
111, refiecting stronger quenched effects. A reverse trend 
is observed for larger bending rigidities. 



F. Discussion of experimental systems 

The following numbers represent typical values for in- 
tegral membrane proteins diffusing in lipid bilayers |27j 
(but are specific to the case of diffusion of band 3 protein 
on the surface of the human red blood cell in the ab- 
sence of osmotic stress and neglecting interactions with 
the cellular cytoskeleton [12]): Dq ~ 0.5 /^m^/s, rj — 0.06 
poise, K/{k^T) ~ 5 and ct = 0. For length scales around 
2'Kq~^ = L ~ 100 nm, thes e n umbers translate to x ~ 3 
and the time scale ratio (57 1 of t^/t'^ — 10^, which 
is close to the annealed limit. However, optical mea- 



surements are limited to (at best) the diffraction limit 
of light, which is closer to 1 fim and cells are typically 
on the order of 10 /im or larger. For these larger mem- 
brane sizes of L ~ 1 — 10 /im, one finds x — 0.3 — 0.03 or 
the time scale ratio of t^/t™ ~ 10 — 1, which spans the 
intermediate regime. While some experimental systems 
may be adequately described by the annealed limit, other 
membrane systems, especially large ones, seem ill suited 
to such a description. The general results presented here 
will be of use in predicting D for systems outside the 
annealed and/or quenched regimes. 



V. CONCLUSION AND DISCUSSION 

We have analyzed the projected self-diffusion co- 
efficient of a free Brownian particle on rough two- 
dimensional surfaces that may be static or fluctuating 
in time. This is achieved by establishing a position 
Langevin equation appropriate for particle diffusion on 
weakly curved membranes represented in Monge gauge. 
The Langevin equation implicitly incorporates curvature 
and gradient effects that hinder the diffusion of the par- 
ticle once projected on a laboratory base plane. We have 
applied this equation numerically (using Brownian Dy- 
namics simulation techniques) to several different exam- 
ples of a static as well as a dynamic surface. 

As shown, the dynamical state of the coupled particle- 
membrane system may be described by a dynamical cou- 
pling parameter, x — / {ArjLDo) (with 77 being the 
medium viscosity, L the system size, and Dq the particle 
curvilinear diffusion coefficient), that controls the sepa- 
ration of characteristic time scales for particle diffusion 
and membrane relaxation. Qualitatively, it represents 
the ratio between the relaxation time of membrane un- 
dulations of typical energy scale k^T and lateral range L 
to the time spent by the particle to diffuse over such a 
length scale. This parameter specifies the dynamical na- 
ture of surface disorder experienced by the diffusing par- 
ticle with two well-defined limiting cases of quenched and 
annealed disorder being given by x ^ (static surface 
or fast diffusion) and x ^ 00 (fast-fluctuating surface or 
slow diffusion), respectively. 

The simulation methods presented in this work enable 
one to study the diffusion process for the whole range 
of non-trivial intermediate couplings Xi where the mem- 
brane and particle diffusion time scales are comparable. 

In all quenched cases considered here, we find that the 
simulation data are best described by the area-scaling 
prediction for the projected diffusion coefficient. Indeed, 
to within our statistical errors, the area-scaling result 
can be considered exact. The area-scaling prediction has 
so far been considered only as a plausible scaling result 
[TBI [5S] and no rigorous proof is available, to our knowl- 
edge, that demonstrates its validity in a systematic fash- 
ion. The variational bounds and the effective medium 
approximation discussed in detail in Ref. |16j show no- 
table deviations from simulation data in all the studied 
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cases. The present findings call for a more detailed anal- 
ysis of diffusion over two-dimensional surfaces, and in 
particular, an analytical explanation that can clarify the 
unanticipated success of area-scaling in the quenched dis- 
order limit. 

In view of the fact that many biological structures ex- 
hibit large static undulations (as in the endoplasmic retic- 
ulum and the lateral cortex of auditory outer hair cells 
EZIHH]), it is important to realize that the annealed re- 
sult may not be applied to these system. Area scaling 
provides the appropriate choice for D/Dq- 

In the annealed limit (x ^ 1), we have shown that 
the simulation data for the projected diffusion coefficient 
agree well with analytical annealed limiting results as 
calculated in Appendix |B] The annealed diffusion coeffi- 
cients are higher than the corresponding quenched values. 
This is more pronounced for strongly roughened mem- 
branes that may be represented by small bending rigid- 
ity (and surface tension) . For vanishing bending rigidity, 
the annealed diffusion coefficient tends to D/Dq — 1/2, 
whereas the corresponding quenched value tends to zero, 
reflecting a huge qualitative difference between the two 
in this limit. For large bending rigidities, the surface 
is almost planar and all different theoretical predictions 
converge to D/Dq = 1, with deviations of only a few 
percent for K/{k^T) > 1. This is why we have exten- 
sively studied the low bending rigidity regime in order 
to emphasize the physical mechanisms which affect the 
diffusion coefficient. 

An important advantage of the present Langevin ap- 
proach is that it can be easily generalized to include 
particle-particle and membrane-particle interactions (as 
in the case of active protein inclusions [IS], where-unlike 
the present case-membrane dynamics are directly influ- 
enced by forces exerted by the proteins). This will en- 
able numerical investigation of more complex and realis- 
tic models of biological membranes. 

As a final note, we comment on some of the physics 
absent from the present study. This work considers the 
influence of surface shape on the diffusion process solely 
through the geometric effect of projection onto a base 
plane. We have neglected any effects whereby surface 
shape may alter the local character of the surface and 
the value of Dq- For example, the curvature of a lipid 
bilayer may influence the value of Z^o itself |25j and we 
have neglected this effect. In the case of fluid lipid bilay- 
ers, it is clear that surface fluctuations must be accom- 
panied by a local flow of lipids as membrane material is 
moved about. Such flows would be expected to drive par- 
ticle motion and lead to correlations between membrane 
fluctuations and protein transport. Our study, and the 
others we are aware of |24l |25], completely neglect this 
effect. The problem of coupled fluctuations and hydro- 
dynamic flow within the bilayer is a seemingly complex 
one, with only formal results available thus far ^49j. 
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Note Added in Proof 

An interesting paper by Reister-Gottfried et al. [50] 
was published while this work was under review. That 
paper also considers the problem of diffusion over an an- 
nealed flexible membrane. 



APPENDIX A: LANGEVIN EQUATION IN 
REDUCED REPRESENTATION 

In this Appendix we derive the rescaled functions that 
enter in the rescaled Langevin equation for particle mo- 
tion, Eq. (51 1. As explained in the text, this follows by 
simply plugging rescaled variables x = x/L, t = Dot/L^ 
and h(x,y,t) = h(x,y,t)/L into the actual Langevin 
equation ( |17[ ) (or simply multiplying both sides of the 
equation by L/Dq). The rescaled drift, Vi, the rescaled 



multiplicative noise factor, fy- (being the square-root of 



the inverse metric tensor, Eq. (21 1), may thus be ob- 
tained as 

Vi[h{{xk},i)] = ^Ui, f,j[h{{xk},i)] ^ Uj, (Al) 

where we have explicitly indicated their dependence on 
space coordinates (particle position) through the surface 
height profile, h{{xk},t) — h{x,y,t) (for fc — 1,2 in 
2D). It is to be noted that the metric tensor gij, and 
therefore its inverse, {g~^)ij, and determinant, g, remain 
invariant under rescaling coordinates, e.g., gij = g. 



Sij + dihdjh, where di = d/dxi. Hence using Eqs. (18 1 



and (21 ) we have 



Sij - 



(g ^)jkdjdkh 
di h di h 



dih 
9 



(A2) 
(A3) 



Likewise, the noise term, rjiit), in Eq. (17 1 is rescaled 
as fii{i) — Lrii{t)/Do. Therefore, fii{i) is a Gaussian 
white noise with zero mean, {fiiit)) = 0, and variance 

A similar scheme can be used to rescale the membrane 
dynamical equation (48) as explained in the text. In 
this case, the rescaled noise term, ^q(f), in Eq. (|52| 
is related to the actual noise term, ^q(i), in Eq. ( |48[ ), 
via £,q{t) = L^q(i) /^/Dq. Therefore, £,q{t) is a Gaussian 
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white noise with {^q{t)) = and variance (^q(t) ^q'(i')) 

2Jq,_q,J(t-?)- 



APPENDIX B: THEORETICAL ANNEALED 
PREDICTION FOR THE RATIO D/Do 

As originally shown by Gustafsson et al. [16 , the dif- 
fusion coefficient ratio D/Dq for a free particle diffusing 
on a fast-fluctuating surface is given by 



D 
Do 



(Bl) 



which has a form similar to the upper bound in the 
quenched limit, Eq. (311, but here brackets denote an 



ensemble average over equilibrium configurations of the 
surface as will be specified below. Gustafsson et al. [T5] 
calculated this average for a generic infinite surface ex- 
hibiting an isotropic Gaussian distribution of surface gra- 
dients, V/i (and thus a scalar diffusion coefficient). For 
an elastic membrane of arbitrary size L (and finite small- 
scale cut-off, a), the gradient distribution is in general 
anisotropic. In this Appendix, we demonstrate how the 
ensemble average in Eq. (Bl I can be calculated for this 
latter case. 

The ensemble average of any function, (5(V/i(p)), of 
surface height gradient, V/i, at a given point p = {x,y), 
can be expressed as a functional integral over equilibrium 
membrane configurations as 



(Q(vMp))),= 



V[h{p')] 



Q{Vh{p)) 



,-/3-H[/i(p')] 



(B2) 

where = l/(fcBT) and Zh = JV[h]exp[-n/ikBT)] is 
the membrane partition function. An example for Q is 
the quantity 1/g appearing in Eq. ( |B1[ ), where g = 1 + 
(Vft-)^. Denoting statistically probable values of V/i by 
u, we obtain a formally equivalent expression for {Q)h at 
a given point p, i.e. 

(g),^ = JduQ{u-p)K{u-p), (B3) 
V[h{p')] 



where 

Puin;p) 



S{u-\7h{p)) 



5(u- V/i(p))e-''^[''(''')l(B4) 



(B5) 



is nothing but the probability distribution function of 
the height gradient, u = Vh, at position p. We thus 
first need to determine Pm(u; p) in order to calculate Eq. 
dBll). 



We proceed by replacing the delta-function in Eq. (B4 ) 
by its Fourier representation, i.e. 



5(u-VMp))=/(^exp 



i/x • (u - Vft(p)) 



(B6) 



where p is a, Fourier conjugate of u. Inserting this into 
Eq. (B4) and using the explicit form of the Hamiltonian 



(37 1 (in Fourier representation), we reach at the following 
expression 

dp 



Pu{u; p) 



(B7) 



with T(p; p) being the so-called generating function of 
height-gradient moments, which is obtained as 

' J 

(B8) 

where I?/iq = Y[q d/iq and fiq = Kq* + aq^ . The above 
Gaussian integral can be evaluated and it follows that the 
dependence on the position of the generating function, X, 
drops out yielding 

(q-M)2- 

q 



I(p- p) = I(p) = exp 



2/3L2 



(B9) 



Putting this expression into Eq. (B7), one finds after 
taking the Gaussian integral over p that 



Pu{n;p) = Puiu) 



g— juA u 



27rVdetX ' 



where the elements of the 2x2 matrix A read 

1 



(BIO) 



(Bll) 



In the limit of large membrane size M = L/a oo, the 
cross terms (off-diagonal elements) of A vanish leading 
to isotropic height-gradient probability distribution used 
in Ref. [TB], i.e. 



P„(u) -> ^ e 



(B12) 



with the variance defined as = ((V/i)^). 

For a finite M = L/a, one can calculate desirable av- 
erages, Eq. (B3l, using the anisotropic height-gradient 
distribution (BIO I. For the particular case of {l/g)h in 
Eq. ( |B1[ ), the computations need to be carried out nu- 
merically. The results for the predicted annealed ratio 
D/Dq, Eq. (|B1[), at different bending rigidities (and zero 



surface tension) are shown in Figures 11 and 12 (dashed 
lines). (Note that the annealed expression for D in gen- 
eral involves additional "off-diagonal" terms containing 
cross-correlation functions, but they are negliggible for 
the membrane size chosen here-see Appendix [C] and Ref. 

Reister et al. [Mj present a method to explicitly eval- 



uate Eq. (Bl I for an elastic membrane, but their scheme 
actually computes {A±/[J d^xg))^^ instead of {l/g)h, 
which in principle represents a different quantity. The 
results reported in Ref. 24J underestimate the annealed 
values of D/Dq by up to 5-10% depending on the bend- 
ing rigidity (with the larger errors being where the dis- 
crete summation over q-modes is approximated by an 
integral) . 
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APPENDIX C: ANNEALED LIMIT: ADIABATIC 
ELIMINATION OF FAST VARIABLES 



(55 1. Inserting Eq. (C4) into the joint Fokker-Planck 



equation (Cl I, one finds 



The coupled membrane-particle dynamical equations 



(51 1 and (52 1 as represented in rescaled units provide an 



appropriate analytical framework for a systematic per- 
turbative analysis around the annealed-disorder limit. 
The time-scale separation in this limit is such that the 
membrane degrees of freedom represent the fast variables, 
whereas the particle degrees of freedom constitute the 
slow variables in the system. Under this condition, one 
can use the methods described in literature to adiabati- 
cally eliminate the fast variables and derive effective dy- 
namical equations for the slow ones. 



Equations (51 1 and (52 1 already indicate how this pro- 



cedure should be done in the present context by making 
use of X as an expansion parameter. Following Risken 
[3|, we can develop a perturbative scheme much in the 
spirit of Born-Oppenheimer approximation. The Fokker- 
Planck equation for the joint probability distribution of 
the set of variables {hq(t), Xi(t)} may be written as 

Pi{S:M,i) ^ {Ox + x^h)Pi{i^,h}^i)^ (Cl) 

where the Fokker-Planck operators associated with x- 
and /i-variables read 



^ dxi ~^ ^ dxidxj 



(C2) 
(C3) 



The idea here is to expand P{{xi,hq},i) in terms of 
the eigen-functions of the operator Oh (which are in the 
present context independent of {ii}), i.e. 

oo 

P{{£,,h^},i) = Crn{{i^},i)(|)rn{{h})■ (C4) 

m— 

Here we have assumed that a stationary solution as well 
as a discrete spectrum of eigen-functions, <t>ni{hq}) (with 
n > 0), exist, i.e. 



Oh (f>,ii{hq}) = -A„ (t)n{{hq}), 



(C5) 



with Ao = being the stationary solution correspond- 
ing to the equilibrium distribution function of membrane, 
Mi^q}) = IlqPcqihq). Notc that Pcq is glvcu by Eq. 



nm ^m: 



(C6) 



m=0 



where Cnm{{xi\) = JVhqcjjilOx (j^m with (f)}^ being 
the eigen- function of the adjoint operator Oj^, which 
is needed since the Fokker-Planck operator is non- 
Hermitian 3 . 

For l arge x 2> 1 , the time derivative may be neglected 
in Eq. (C6 1 for n > 1 in the long-time limit t ^ (xAi)~^, 
and only kept in the equation governing cq (with Aq = 0) . 
By doing so, one can easily derive a perturbative expres- 
sion for the latter equation. But it turns out that Cq is 
nothing but the time-dependent probability distribution 
of the slow variable averaged over the fast-variable dis- 
tribution [3], i.e., Co = 'P{{xi},i) — J'DhqP{{xt,hq},i). 
Therefore, one ends up with the desired effective equa- 
tion governing the probability distribution of particle po- 
sition on a fast-fluctuating membrane, which admits the 
following perturbative form 

viix.}, t) = (^61") + oi'^ + o{x-'?j V{{x.}, t), 

(C7) 

with the zeroth-order term being given by 



6(°) ^ Cooi{S:^}) 



_d_ 
dxi 



{Vi)h ■ 



dxidxj 

(C8) 

inverse metric ten- 
~^)ij)h, are averaged 



where the drift, {vi)h, and the 
sor (rescaled diffusion matrix), ((g 
over the equilibrium distribution function of membrane 
height at a given point as, e.g., d efine d in Eq. (B2| 

The first-order operator in Eq 

1]^=! Gondii}) Cno{{Xi})/Xn. 



( C7 1 reads 



Note that the zeroth-order term in Eq. (C7) becomes 
exact in the annealed limit x — > cxd. Using explicit ex- 
pressions given before for Vi and (g^^)ij (and neglecting 
small off-diagonal terms for large L), the annealed ex- 
pression for D/Do, Eq. (Bll, immediately follows. The 



effective Fokker-Planck equation for particle diffusion on 
an annealed surface was first derived in Ref. [24 using 
a pre-averaging approximation, which gives the effective 
Fokker-Planck equation ( C7 1 up to the zeroth-order term. 



The preceding derivation places this result within a sys- 
tematic framework and predicts higher-order corrections. 
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